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An explicit representation of the transition densities of the skew 
Brownian motion with drift and two semipermeable barriers 
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Abstract 

In this paper, we obtain an explicit representation of the transition density of the one-dimensional skew 
Brownian motion with (a constant drift and) two semipermeable barriers. Moreover we propose a rejection 
method to simulate this density in an exact way. 
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1 Introduction 

The need to study the skew Brownian motion, and in particular its explicit transition densities, has emerged in 
various contexts during the last years. An overview and list of historical background and main applications can 
be found in [1]. Nevertheless, to the best of our knowledge, the transition densities of the one-dimensional skew 
Brownian motion with constant drift and two semipermeable barriers was not yet given as a closed formula, not 
even for the driftless version. In the latter case, one can only find in [S] a non explicit formula for it. 

We obtain here a closed formula for the transition density of the skew Brownian motion with drift and two 
semipermeable barriers as series of Gaussian transition densities and cumulative distribution functions. This is 
a non trivial generalization of the case of reflecting barriers treated in |2Ij . 

In order to avoid repetitions, from now on we will use the following notation: /3-SBM is the skew Brownian 
motion with one semipermeable barrier of permeability coefficient /3 and (/3i,/32)-SBM is the skew Brownian 
motion with two semipermeable barriers of permeability coefficients respectively /3i and ^ 2 - 

The /?—SBM was introduced by Ito and McKean in [10], as a one-dimensional Wiener process transformed by 
flipping the excursions from the origin with probability G (0,1) (if /? = 0 it is the usual Brownian motion). 
Unfortunately this trajectorial definition does not lend itself to generalizations. 

The skew Brownian motion behaves as a Brownian motion between the barriers but it has a particular 
behaviour when it reaches them: it is partially reflected. This interpretation yields the various generalizations, 
that we are going to present shortly. 

A recent survey on the skew Brownian motion can be found in m in which various equivalent representations 
of the semigroup are given. Let us now present the process as a solution to a stochastic differential equation. 

It was proved by Harrison and Shepp in |5| that if |/3| < 1, there is a unique strong solution to the stochastic 
differential equation involving the symmetric local time at the point 0 {L^)t>o 

idXt=dWt+/3dL°{X), 

\Ao = a:, L? =/A{v,=o}dL°, 
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that is the /3—SBM. In particular if /3 = 0, it is the usual Brownian motion. Harrison and Shepp also proved 
that if |/3| > 1 there is no solution to 0. Notice that if a: > 0, the 1—SBM is the reflected Brownian motion on 
the positive semi-axis, and if a: < 0, the (—1)—SBM is the reflected Brownian motion on the negative semi-axis. 


There are many possible generalizations of the SBM: one-dimensional skew BM with more semipermeable 
barriers m ,[14],[18]), n-dimensional skew BM with one permeable barrier, as it is called in [14] referring to 
m and m, and distorted Brownian motion m)- A new proof of the weak existence and uniqueness for the 
n—dimensional SBM appeared recently in [5]. 


The existence of several barriers does not allow anymore a trajectorial interpretation as randomly flipped 
excursions like for the ,5-SBM, nevertheless one can define the process as the unique strong solution to a slight 
modification of equation 0. The stochastic differential equation E ((,5i,/32),/r) satisfied by the (,5i,/32)—SBM 
with drift p, G K is indeed 


dXt = dWt -b /rdf -b /SidTf (X) -b {X), 

Xo = X, Ll^ = 


^((/3i,/32),m) 


where the coefficients /3i,/32 S (—1,1) and Z\,Z 2 G K are the barriers. Obviously, if /32 = 0 the second barrier 
disappears and one obtains the equation satisfied by the /3i-SBM with drift with semipermeable barrier zi: 


dXt = dWt + ndt + PidLt" (A), 

Xo = x, =/ol{x.=zi}rfir. 




The transition probability density function x, j/) of the Markov process, unique solution to f(/3i,^) 

is computed in [^ using the trajectorial interpretation. As already noticed this approach is not extendable for 
finding the transition density in presence of more barriers. So let us briefly recall how to compute the semigroup 
of the /?—SBM with barrier in zero as solution of a partial differential equation with specific boundary conditions. 
In [16] and m it is shown that 

( 2 ) 


L — - A -b /3doV 


is, formally, the infinitesimal generator of the /3— SBM with barrier in zero. Moreover the parabolic problem 
dtu = Lu (whose solution is the semigroup generated by L) is equivalent to the transmission problem (see [H], 
section 3.1): 

f dtv = \Xv 

I (1 -b (3)Xv{t, 0+) = (1 — P)Xv{t, 0“) (transmission condition). 

A solution to (|^ is equivalently a weak solution to the following problem: 


( 3 ) 


f uit, x) G C(0, T; L2(K)) n L2(o, T; ^^(K)), 

< dtU = Lu, (4) 

[m(0,x) = v?(x) G T^(K), 


where L is the divergence form operator 

with domain I?(L) = {(b’S k{x)(p'{x) G iJ^(IR)}. Using Dirichlet forms one proves that the unique 

solution of 0 is the semigroup of the /?—SBM, solution of 0 (see section 3 of [H]). 

Our approach for computing the transition density of the (/?i, /? 2 )—SBM with or without drift will be based on 
identiying its infinitesimal generator as a divergence form operator, generalising the case of the driftless process 
treated for example in [S].|13j. Once we will have computed the divergence form of the infinitesimal generator 
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{L,V{L)) associated to g((/3i,/32 ), m)I solve the Kolmogorov equation satisfied by the semigroup: for 

each continuous and bounded function /, P^f is the solution in x M \ {0},M) nC(M_|_ x IR,IR) of 


§^u{t,x) = Lu{t,x) = ^^u{t,x) + t e K+, a; e K\ {zi,Z 2 }, 

l+P: 


2 Vu(t, Z+ ) = Z~ ) 

L{t,Z+) = u{t,Z~) 
i{0,x) = f{x) 


t e ] 
t e ] 

X G 


^+; J — 1) 2, 

^+; j — 1) 2, 


(5) 


The transition density {t,y) i-G p{t,x,y) will satisfy, for x fixed, the analogous PDE for the adjoint L*: 


§iPit,y) = L*p{t,y) 

\\'p{t,z+) - ppit,z+) = \Vp{t,z-) - pp{t,zj) 
(1 + I3j)p{t, zj ) = (1 - Pj)p{t, z +) 

,p(0,?/) = 4(2/) 


t G (0, +oo), 2 / G M \ {zi,Z 2 }, 
tGM+, j = l,2, 
teR+, j = l,2, 

2/ G K. 


( 6 ) 


The paper is organized as follows : in Section we give an explicit characterization of the infinitesimal 
generator associated to the solution to £^((/3i, /32), y) in order to obtain a representation of its transition density. 
Then we exploit it in the following cases: first for the (/3i,/32)—SBM without drift, then for the /3—SBM with 
constant drift, and finally we give the formula for the drifted version of the (/3i, /32)—SBM. Moreover we discuss 
some particular and limit cases and compare our results to former ones. In Section we present a rejection 
sampling method that allows to simulate exactly the SBM with two semipermeable barriers. 


2 The transition density of the (/3i,/32)-SBM with and without drift 

2.1 The framework and the method 

In order to obtain the transition density of the (/?i,/32)-SBM, we identify its infinitesimal generator. The 
infinitesimal generator of the /3i-SBM with one semipermeable barrier in zi, solution of the equation f(/3i,0), 
is the divergence form operator 

{ml) , m) = {V' e KxW{x) g h\r)] , 

k{x) = i + 4 (l[,,^_+oo)(a;) - \) 

with piecewise constant function k{x) unique up to a multiplicative constant. (See for example [13jb 

Notice that a straightforward generalization of Q yields to the generator of the /? = (/3i, 4,. ■., ,0n)-SBM 
with n semipermeable barriers in zi < Z 2 < ... < Zn, by modifying the piecewise constant function k{x). Indeed, 
in case of two semipermeable barriers zi,Z 2 j the function k{x) assumes three different values: 

\{l - j3i){l - (32) x<zi, 

3(1 +/ 3 i)(l - ^2) zi<x<z 2 , (8) 

3(1 +/ 3 l)(l + ^2) X>Z 2 - 

Therefore taking k = nm=0 (2 (l[zm,+oo) 2 )), {L,'D{L)) is the infinitesimal generator of the SBM with 

n semipermeable barriers. 


m = +/^i (lbi,+oo)(a:) - 


. I 


The operator L is a divergence form operator with discontinuous coefficients, therefore one obtains a repre¬ 
sentation for the transition densities, as in [5], Chapter II. In section 5 the authors recover themselves the case 
of the /3-SBM and the (/3i,/32)-SBM without drift. Unfortunately in the latter case the authors do not explicit 
further the transition density function, they just identify it through a “kind of 0-function” 

1 r +00 

h{t,i,C,a) = — 6-“ (1-b 176™“) du>. (9) 

27r 7-00 
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We are therefore going to generalize the method, called Green function method or Titchmarsh-Kodaira-Yoshida 
method, for giving an explicit representation of the transition density function associated to the slightly more 
general infinitesimal generator including a constant drift ^ S K. 


I?(L) = {'i/i S Hl{h{x)dx)', h{x)'tjj'{x) € ff^(k~^(x)dx)} 
^h(x) := k{x)e‘^f^^ 


( 10 ) 


where k{x) is the piecewise constant function defined in 
bounded from above. 


Remark that h{x) is strictly positive but not 


Lemma 2.1. The operator (L,I?(L)) defined by (10) 

(i) is self-adjoint in Lf(h{x)dx) and its spectrum a{L) is a closed subset o/(—oo,0] containing 0; 
(ii) is the infinitesimal generator in L‘^{h{x)dx) of the process [j5i, j32)-SBM with drift /i. 

Proof. First of all notice that the measure vi^dx) := h{x)dx is not a finite measure. The form 

q : Hl{h{x)dx) x Hl{h{x)dx) —)■ K defined by (/, 5 ) i-J- f f'g'h{x)dx, 

Js. 


( 11 ) 


is symmetric, semibounded and closed with domain Q{q) = Hl{h{x)dx) C L'^(h{x)dx). Therefore there exists a 
unique operator T with 'D{T) C Q{q) such that q{u,v) = — < u,Tv >v (see for example Corollary 1.3.1 in [3). 
Moreover the operator T is self-adjoint. 

One can show that {2L,'D{L)) = {T,'D{T)) using the implicit characterization of I?(T). Therefore the operator 
2L is self-adjoint, hence the conclusions on its spectrum. 

We now apply the results presented in the recent paper |15j , Remark 2.6.ii: the Hunt process whose semigroup 
is associated to the closed form {q, Q{q)) in 0 is the SBM with drift with semipermeable barriers. By 
uniqueness of the self-adjoint operator associated to the form we conclude that the operator {L,V{L)) is the 
infinitesimal generator of this SBM. □ 


Remark 2.2. (i) The same lemma holds for the fli-SBM with drift, and also for the driftless processes 

{Pi, h)-SBM and Pi-SBM (p = 0); 


(ii) As an alternative to one can express the infinitesimal generator for the {Pi, P 2 )-SBM with drift as 


V{A) := {if e H^{dx)- kip e H^{dx)} 


with k{x) defined in Q). In that case, A is not self-adjoint. 

(Hi) One can show, using Hille-Yoshida theorem, that the operator {L,'D{L)) is sectorial since it is self adjoint 
and in particular it is the infinitesimal generator af a strongly-continuous semigroup of contractions. 


Since the infinitesimal generator {L,T>{L)) is a sectorial operator, its associated transition semigroup P* can 
be represented as: 

PMx) = ^ J e^*ux^,f,{x)dX 

where T is a contour in the complex A plane around the negative semi-axis (— oo,0] that contains the spectrum 
cr{L) and is the resolvent solution to (A — L)u\^,p = p for all Lp G L‘^{h{x)dx) (see for example Theorem 
12.31 in [H]). Therefore the transition density satisfies 


P{t,x,y) 


1 

2Tri 


e^*G{x, y, X)dX 


( 12 ) 


where G{x, y; A) are the Green functions. 

^The authors would like to thank Markus Klein (Universitat Potsdam) for the interesting discussions 
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Lemma 2.3. For each A S C \ K_ the Green functions are given by 


G{x,y]X) 


-2h{y) 


U-\-i^y , X')U— (x, A)l|a;<y| -j- [/_!_(x, A) U— (y, 

/i(xo)PL(J7_,[/+)(xo,A) 


where U± S 2A(L) are the solutions to: 


(A-L)t7±(x,A) =0, 


lim C/+(x, A) = 0, 

X—>- + 00 


lim U- (x, A) = 0, 


(13) 


while W{U-,U+){xo, X) = {7_(xo, A)17(_(xo, A) — U'_{xo, X)U+(xo, X) is the wronskian in xq S K. 

Proof. One can easily prove that the function x i-)- h{x)W{U-,U+){x, X) is constant and check that x i-)- 
G{x,y]X) is a solution to (A — L)v{x) = 6o{x — y) for all y G M. By uniqueness of the solution the proof is 
done. □ 

Remark 2.4. Notice that x i-G G(x,y,X) G T){L), and y >->■ G{x,y;X) G V{L*) = G I?(L)| since 

L*g := hL (^) is the adjoint in L^{dx). The same holds for y i—>■ p{t,x,y) G T>{L*). 


2.2 The case of (/3i,/? 2 )-SBM without drift 

We will now present the method step by step. 


2.2.1 The Green functions 


The first step is to find the eigenfunctions 17+(x, A) and 17_(x, A) of L defined in (13). The two barriers divide 
the real line into three intervals over which the functions U± can be constructed as linear combination of the 

eigenfunctions of the operator L for the eigenvalue A G C \ (—oo,0] that are m_(x) = exp M+(a;) = 

exp v^x'). Therefore 


and 


U- X < Zi, 

U- = A{X)u- + i?(A)u+ 2:1 < X < Z 2 , 

C{X)u-+ D{X)u.^- X > Z 2 ] 

G(A)m_ + 77(A)m+ X < zi, 

U+ = ^ E(X)u- + F(X)u+ zi < X < Z 2 , 

M+ X > Z2, 


with eight coefficients to be determined. Notice that since U± G 2A(L), they are continuous functions and have to 
satisfy the so-called transmission conditions derived from the continuity of x 1 —>■ k{x)U±{x, A). These conditions 
will determine uniquely the eight coefficients: 


^(A) — (1 -I- /3i) 

B{X) = A{X)l3ie^^^^ 

' C(A) = + 1) ((1 + ^i)(l + /32))”\ 

77(A) = + /32e2^"=) ((1 + ^i)(l + P2)r' 

'g(A) = - (/32e-2^^= + ((1 - ^,){1 - 

77(A) = (/3i/32e-2^(^=-^i) + l) ((1 - ^i)(l - ^ 2 ))”', 
E{X) = -F(A)/32e-2^^G 
F(A) = (1-/32)-i. 
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The second step is to compute the wronskian. Consider xq < zi, for example xq = zi — 1. Hence the 
wronskian is 


WM = = -2v^H(A). 

where we will denote by 0 the distance between the barriers Z 2 — Zi. This leads to the following 
Proposition 2.5. The Green functions are given by 


G{x,y,X) 


where (j>{X) := z := Z 2 — is the distance between the barriers, and 


f ci(?/,/3i,/32) = 1 

I C2(2/,/3i,/ 32) = (21[zi,+oo)(y) - 1) Pi 

I C3{y,Pi,P2) = (21[22_+oo)(y) - 1) P 2 
[c 4 {y,Pl,P 2 ) = (1 - 21 [ zi , z 2 )( j /)) P1P2 


ai{x,y) = 0 

a 2 ix,y) = \y - zi \ + \x - zi \ - \y - x\ 
a3ix,y) = \y- Z 2 \ + \y- Z 2 \ -\y-x\ 

Oiix, y) = 2 {z 2 - max{x, y, zi))’*' + 2 {min{x, y, Z 2 ) - zi)’*' 


Proof. We will only do the computations in the case x < zi < Z 2 < y, the other cases are similar. From Lemma 
|2.3[ since h = k and chosing xq < Zi, the Green function is of the following form 


G{x,y;X) 


-2k{y) 


U+{y,X)U_{x,X) 

k{xo)W{U-,U+){xo) 




1 4:k{y)e 'Xxx{y x) 


It is then sufficient to check that aj{x, y) = 0 for j G {1, 2,3,4} and P) = 4k(y). □ 

Remark 2.6. fij The function p is well defined as bijection between C \ (—oo,0] and G C; 3?(C) > 0}. 

(a) The denominator A i-G 1 + PiP 2 e~'^^^^">^ has no zero m C \ (—co, 0] since ^(j){X) > 0. 

(Hi) aj{x,y) >0 for j G {1,2,3, A}. 


2.2.2 The transition density as (contour) integral 

Since the Green functions depend on A only through 4>{X) = we can apply the change of variable 

A !->■ (j){X) =: ( to the integral appearing in (12): 




e 2 


G{x,y;p{X)) #(A) = 


/ 0 (r) 


e 2 ‘ G{x,y;0 


2.5 


where G{x,y; ({/{X)) = (p{X)G(x,y, X) and G{x,y,X) given in Proposition 
^2 _ 

Since the integrand e^*G{x, y; is holomorphic in the closed subset of the complex plane between iK and 4>{T), 
we could deform (shrink) the contour ^i>(r) to the imaginary line by an homotopy. Indeed, if we denote by M 
the unique point with imaginary part u in ^(T) (as in Figure]^, it is possible to shrink the contour (/)(r) to 
if the following lemma holds: 

Lemma 2.7. Consider the function 


G(x V £) - 

/3i/32e-2«^ + 1 


Then 


lim 

U—>-±oo 


e 2 ‘ G{x,y;(,) d( = 0, 


pM 


where pm is the segment in the figure connecting the point M with its projection on zIR, M' = (0,u). 
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Proof. Let us show that the absolute value converges to zero: 


e 2 ‘ G{x,y;f) df 


< 


e 2 * G{x,y;^) 


df = 




e 2 


G{x,y;iu + v)\ dv 


with i{u) := \M' — M\ (hence M = {£{u),u)) and limu_>oo ^(u) = 0. Let us notice that 


\G{x,y,iu + v)\ < e —- 


|;5i/32e-2(“+2^)2= + l| 


1- |/3i/32|e-2-- 


therefore 


e 2 * G{x,y;f) df 




< e 2 - / 6 2-6 

Jo 




dv 


that clearly converges to zero if |u| goes to infinity. 


□ 




Figure 1: 

a) The picture shows the green image of the blue contour F under </). The spectrum of the operator [L,'D[L)) 
is contained in the red semi-axis (—oo,0], which coincides with the complement of the domain of </>. 

b) In this figure one sees the magenta segment pm connecting the unique point M in (^(F) with imaginary part 

u to its projection M' on the imaginary line. The homotopy iL : [0,1] x K —)■ that deforms (/l'(r) into iK 

is given by H{t, u) = M'(l - t) -f tM. 


Therefore the integral in (12) becomes (with ^ = iw) 


pv 




/3i/32e 


— 2iwz 


dw. 


(14) 


One can also rewrite it using the functions h defined by equation 

= ^Cj{y,l3i,P2)h {aj{x,y) + |a; - y]),/3i/32,, 

j=i V / 


which agrees with the results in [8]. Nevertheless, since |/3i/32| < 1, we can explicit further the expression (14). 
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2.2.3 The transition density as series of Fourier transforms 

Proposition 2.8. The transition density of the {/3i, /32)-SBM has the following expansion 


p^^^'^'^\t,x,y) =p^°'°\t,x,y)^{-(3i(32f^Cj{y,Pi,l32)e ^' 2 * ^ g * (15) 

fc=0 3=1 

where {t, x, y) is the transition density function of the Brownian motion. 


Proof. Let us consider the expression (14). The denominator can be seen as the sum of a geometric series 


1 

1 + /3i/32e-2™^ 


00 


since |/3i/32e = \PiP 2 \ < 1- 

Therefore the density can be written as 


1 2 

^ fc=o i=i 


We can exchange integral and series, because the series of absolute values e~~* is integrable. 

We conclude that the transition density is a series of Fourier transforms: 


°° ^ 1 /■ 2 

= V(-/3i/ 32)" Vc,(y,/3i,/?2)— / 

fc=0 i = i 

- CXD 4 

= +2zA: + |x-?/|) 


(16) 


/c=0 

00 


i-1 

4 


1 V^/ n n \kSr^ / o /D ^ + lx - 2/|\ 


y/riri 

where gtiw) := = gi{w\/i) and its Fourier transform satishes gt(w) = ^gtij) = ^9i (^) • 

We notice that 51(0 + b) = gi{a)gi{b)e~°'^ hence we can write the density as 




\/27rt 


fc=0 J—1 


Using the identity x, y) = conclude and obtain (15). 


□ 


2.3 The case of the (/3i,/32)-SBM with drift 

2.3.1 Expansion of the transition density in the case of one barrier and drift 

In this subsection we propose an explicit computation of the transition density function of the /3i-SBM with 


constant drift /r, solution to £ (/3i, p) 


Proposition 2.9. The transition density for the SBM with constant drift p and barrier in zi satisfies 


X, y) = p''^^ {t, X, y)v\f^\t, x, y), 


r.(o) 


{Mi 


















where p^^\t,x,y) is the transition density of the Brownian motion with drift y (without skew), and 

:= U - l{xi!,i>o} + [l + A (21[;i,+cc)(p) - l)] exp f-V 

^ 1 - ftpvstexp + lai+ , 

OO 

where xi := x — zi, yi := y — Zi and := e~^du is the queue of a standard Gaussian law. 

This result appears in [5] for a barrier in zero although it holds for any barrier. The authors prove it using 
the trajectorial definition of the SBM. 

We provide here a completely different proof based on the generalization of the Green function method. Indeed 
the infinitesimal generator of the process is a generalization of Q with the function h{x) = k{x)e^^^ instead of 
k{x). We will denote by (3 := /3i G (—1, f) \ {0} the unique skewness parameter. The same method we develop 
here, will also provide the transition density for the (/3i, /32)-SBM with drift even though trickier technical issues 
are involved. 


A. The Green functions 


When there is a drift y 0, the functions U±{x,X) solutions to (131 are linear combinations of the two 
eigenfunctions for the eigenvalue A S C \ (—oo,0] u±{x,X) = exp (^—yx^f y^y'^ + 2Xx'^. The coefficients are 
uniquely determined using the continuity and the transmission conditions: 


u- on(-oo,zi], ^ _ I £;(A)m_ + F(A)u+ 

A{X)u-+ B{X)u+ on[zi,+oo), lu+ 


on (- 00 , zi], 
on [zi, +oo). 


with 

IAl(A) = ^ + ^7^^) > BiX) = (1 - A(A))eV^^-i 

^ (^1 + . E{X) = (1 - F(A))e-2y^-i. 

We compute the wronskian at the point Xg < zi and obtain 

W{U-, U+){xg, X) = —2y/ y? + 2AF(A) exp {—2yxg). 


This leads to the result: 

Lemma 2.10. The Green functions satisfy 

U_{x Xy,X)U+{xy y,X) 


G{x,y-X) = 2h{y)- 


= guiy-x) 


Py + yj 2X + yf 

\ g--\/2A+/i,2|y_a;| 




yf ^2X + yf 


Py 


^Cj{y,y, ^2A + /r2)e-y^AT^a,(x.y) 


U'=i 


where 

\ci{y,y]w) = py + w 

\c 2 {y,y,w) = Pw (21[^j_+oo)(y) - l) - Pt, 

Notice that 02 ( 2 ;, y) > 0 for all x,y gR. 


ai{x,y) = 0 

02 ( 2 :,y) = |y-zi| + | 2 ;-zi| - |y- 2 ;|. 


(17) 
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B. The transition density as a contour integral 

The dependence on A of the Green functions given by ([T7|) is actually a dependence on 


0(A) := x/ 2A+7? g e C; 3fi(C) > 0}. 


This allows the change of variables ^ := 0(A) as in the subsection 2.2.2 

pi^K't,x,y) = X)dX = 




-y/2\+y^\v-x\ f 2 


27rj Jt a/2A + + 0^ 1 


x/2A+7I^)e-V2A00^a,(x.y) j 


dA 


\/2A +, 




2 g-C|y-a:| / ^ ^ , 

0 1 - I dt 

27r^70(^) ^ + 


0 /x > 0 





b 2 ) 


Figure 2: 

a) The picture shows the green image of the blue contour F under 0. The red line (—oo,0] contains the 

2 

spectrum of the operator (L,I?(L)). The dashed line (—c», is the complement of the domain of 0. 

bl) Case /3/x > 0: the figure represents the magenta segment pM connecting the unique point M in 0(F) with 
imaginary part u to its projection on the imaginary line. The red segment (0, |^|] is the image under 0 of 

b2) Case /3/r < 0: the curve 0(r) is decomposed as the union of a green curve 0' that avoids the unique pole 
—fdp G (0, |/x|] and the blue cycle 7 containing it. The segment pM connects in this case the unique point 
M in 0' with imaginary part u to its projection on the imaginary line. 


If fdp > 0 the integrand is holomorphic on the region between the contour 0(r) and the imaginary line. If 
I3p <0 the integrand has exactly one pole of order one in ^ = —/3p. We then decompose the curve 0(r) as the 
union of a curve 0 ' and 7 , where 7 is a loop around the pole and 0 ' avoids the pole. 
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Respectively (/>(r) and </>' can be deformed to the imaginary line (through H : [0,1] x M —^ given by 

H{t,u) = M'{\ — t) + tM), if the analogous of Lemma 2.7 is satisfied: 


Lemma 2.11. 


lim 

|ii |—>+00 , 




i=i 


IpM ? + 

where pm is the segment in the figures connecting the points M' and M. 
Proof. First of all notice that this integral is equal to 


lu := 


e-i\v-^\ I J. 


' pM 




] df = j e-f"-"' ( 1 + 




' pM 


C2{p,y;0 


y) 


df, 


with £(u) := \M — M'\ and lim„_v„£(u) = 0. Let us consider the following parametrization of the segment 
p{M) = {M' + t;(l, 0); vG (0, ((u))}, then: 


\lu\ < 








C2{p,y\v + iu) 


V + iu + /3/i 


^-va2{x,y) \ 


For u large enough 


therefore 


C 2 {p,y]v + iu) 


V + iu + Pp 


I / (^( 21 [zi.+oo)(y)-l) 

'V {v + pp)^+u^ - ’ 


141 < J e - (1 + dv, 

that converges to zero if |m| goes to infinity. 

We compute the integral on the loop through the method of residues: 


□ 




1 f 


Y,c,{p,y-Oe-^^P^’y'i df 


27rf 4 £.+ Pp 

1 


0 ' ^ + Pp \ 




(18) 




E 

i=i 


1 


27r Jr iw + Pp 


e-^‘ 


Cjip, y; 


-Pp{l + P {21[,,,+^fy) - 1)) 


(♦) 


The last equality is obtained by shrinking J^, —> Jfjj and changing variable f = iw. 

C. The transition density as a sum of Fourier transforms 


We interpret each of the two integrals in the last equality of equation (181 as the Fourier transform computed 
at the value {aj{x, y) + |a; — y|) of the function 


w H> r- -^Cj{p,y;iw) = 


if j = 1 


_C ' ( // U' 'I'lJJ ) \ / \ 2 

iw + Pp^^^^^^ |^(^(2l[,^ +^2(y)-l)+;,(l + /3 (21[,^,+^)(y)-l))2^je-^‘ ifj=2, 
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where the Fourier transform of / is J^{f){uj) = f{uj) = /^e “^/(y) dy. In both cases, these functions are 

integrable in ic, so the transition density can now be written as 








E-^ 


i=i 


1 


iw + Pfi 


Cj{y.,y;iw) {aj{x,y) + |a: - y\) 


(19) 


- /3m (1 + ^ (21[.„+oo)(2/) - 1)) 


We can assume that /3/i 7^ 0 because if /3 = 0 we get the simple Brownian motion with drift without skew, and 
if /r = 0 we get the /3-SBM whose transition density is already known (see for example m)- 


Lemma 2.12. If a £ 


then 

T 


(w) = iV^ (2lR+(a) — 1) 1 r- (auj). 


Proof. It is true since 

1 


= J- 


^ (iy/^ (2lR+(a) — 1) 6““^ 1r- (aw)') (w) = — E [ (iV^ (2lR+(a) — 1) e°‘‘^ 1r- (aw)') 6*“^“ dw. 

' y/2'K Jr V ^ 


Notice that is not integrable but iy/^ (2lR+(a) — 1) e“‘^ 1r- (aw) is integrable. 
Using F (‘^) = Lemma 

^tC2{ti,y]iw)' 


□ 


2.12 


we get 


F e"^ 


iw + PyL 


(w) = 


(21pi.+oo)(y) - 1) e - ^1/3^1 (l + /3(21[^^_+<^)(?/) - 1)) • (^/rw) *e ) (w). 


yft 


We compute the convolution as 




{.Ppw)*e ) (w) = v^e''' 2 ’ j %_ (^^) + (21r+ (/3^) - 1) ^ + /3/rv^ ) j. 


I S. 


(* 


Notice that the term (**) arising from the convolution is actually opposite to the term (*) in (18) arising 
from the integration on the cycle 7 containing the pole. Therefore the transition density becomes 


P^i^\t,x,y) = '* 2*^6 +/3(21[^^^+^)(y) - 1 ) e” 


(a 2 (x.y,zi) + |x-y|)^ 


- (1 + /3 (21[^^,+oo)(?/) - 1)) 


a 2 {.x,y,zi) + \x-y\ 

Vt 


Ppy/i 


Isolating the density of the Brownian motion with drift y, without skew, we recognise the expression we wanted. 


which completes the proof of Proposition 2.9 


(\x-y\£-\x-y£ 


plfHt,x,y) =P^°Ht,x,y) [e +/3(21[^^_+^)(y) - 1 ) e’ 


(\x-zi\ + \y-zi\)^ -Ix-yl-^ 


+ (2I|..,+„,(„) - (^EzjA±i]L_F±M '^. 
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2.3.2 The case of two barriers and drift 

In this subsection we extend the computations done in the previous one to the case of two barriers to provide 
the transition density for the (/3i,/32)-SBM with drift. 

Theorem 2.13. Suppose Pipt > 0 and P 2 P- > 0. The transition density of the {Pi, P 2 )-SBM with drift decomposes 
as 

X, y) = (t, X, y)vf^’^^\t, X, y) 

where the function is given by a series of Fourier transforms. If Pi P 2 , 


,(/3i./32) 


{-P 1 P 2 )'" ^ ^ (-l)’"(2fc-n)! ^ 


],2-h{y) 

k=0 ^g^(j(fc-n)!(fc-m)!n!m!(/3i-/32)2fc-"2^^ {nyp)h 


EE Trli; 


_V V 

k—n / ^^ 


and if pi = P 2 , 


( 20 ) 


,(/3i./3i)/ 


{t,x,y) = e^ Y. EE E ( 

( 21 ) 


4 2 fc 


{2k + 1)! ^^—' \m 

k=0 ' ^ j = lh=0m=0 


where ojj^k ■= z ■.= Z 2 — zi and aj{x,y) and Cj^h{y) o-re defined in Lemma 


2.14 


:=<i^Y(uj,P2pVi) - (-l)”^^_„(a;,/3i/rVt) 


and for A G {Pifiy/i, P 2 pVi}, 


'"+LtJ (i)i+h 

\Juj,A) = {2m + h)l Y M 


e=o 


2^ t!(2(TO-£) + h)! 


where 


and 


n 2(m ^) + /i . . /o/ I I.\ 

SY.i(‘^,A)=Y E (j( )(u; + 24)"-M2(—^)+'*-V,+.(a2,^), 

IT*—n a —n \ X \ / 


r-=0 s=0 


'^e^+^‘^<l>^{uj + A) 


Jq{uj,A) := < 


—e 2 


9 = 0 , 

9=1, 


Jo{uj,A){q-1)!!-Ji{u;,A)J2^^g{u; + A)'i q>2 even, 


9-1 v 


r-^-o 


q > 3 odd . 


The proof of the theorem is based on the following four lemmas. 
Lemma 2.14. The Green functions, defined in Lemma\2.3\ satisfy 


(Jtq. ,,.,,.'1 _ Ay-x)p-w\x-y\ _ Z^j_i 3 \P^y^ ) _ 

^yx,y,w, ^e -/r2 ) + (^ + ^^^)(^ +/j^m) ' 


where w := i/2X + p'^ and z := Z 2 — Zi is the distance between the barriers. The functions aj{x,y) are non 
negative, in particular they are 

' ai{x,y) = 0 

a- 2 {x,y) = \y- zfi + |a: - zi\ - \y - x\ 
a-z{x,y) = \y- Z 2 \ + \y- Z 2 I -\y-x\ 

^a4{x,y) = 2 {z 2- max{x,y,ziY + 2 {min{x,y,Z 2 ) - zY 
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and Cj{fj.,y,w) = w'^Cjfi{y) + wy,Cj^i{y) + y?Cj^2{y) where 


ci,o(y) = 1, 

C2,o( 2/) = (21[zi,+oo)(y) - l) Pi 
C3,o{y) = ( 21 [ z 2 .+ oo )( 2 /) - 1 ) P 2 
Si,o{y) = (1 - 21 [ zi . z 2 )( 2 /)) P 1 P 2 
Proof. Analogously of the proof provided in subsection |2.2.1| 
Lemma 2.15 (Partial fractional decomposition). Let a,b Gl 


ci,iiy) = P1+P2 
C2,i{y) = -Pi - C4,oiy) 
C3,liy) = -P2 + c4,o{y) 
C4,i{y) = 0 


'ci,2{y) = P1P2 
C2,2{y) = PlC 3 ,o{y) 
03,2(2/) = -P2C2fi{y) 
.04.2(2/) = -04.0(2/)- 


□ 


a b, then 


{w — za)^+i(w — ib)^~^^ 


k 


i=o 


(a — 6)2fc+i-i 


2 k - j 
k- j 


{w — z&)-^+i 


- (-1)^ 


{w — m)f+i 


Proof. The function f{x) = is a rational function with two poles Xi = ia,X 2 = ib of order 

fc + 1. We followed a standard method for computing the decomposition: there exist coefficients 0;^^ such that 
the function can be written as f{x) = X)i=i {x°^x y ■ ^iace the ttij are the residues in Xi of the function 

gi,j{x) = (x — Xiy~^f{x), we computed them explicitly. □ 


Lemma 2.16. If a 

T 


and fc € N then 
1 


{w — m)''+i 


^ (w) = (2lR+(a) - 1) ^ e°-‘^ 1r- (aw). 


Proof. If fc = 0 it coincides with Lemma 
computes its Fourier transform in w, fj 


2.12 


otherwise the function 


e Li 


n 


(w — ia)^+^ 


"dw, through the method of residues. 


Lemma 2.17. Let q £ N. The primitive function Iq{-) of v ^ v^e~^ (resp. Iq{a) = 
satisfies 

Io{a) := •\/^$(a) = {resp. •\/27r$'^(a)) 

^2 ^2 
/i(a) := — (resp. e~^) 

Iq{a) = Ioia)iq - 1)!! + Ii{a) 


and one 

□ 

^2 

e~^dv ) 


Iq{a) = 


q = 0, 

9 = 1, 
q > 2 even, 
q > 3 odd . 


Proof. Straightforward for q = 0,q = 1, and for q > 2 one can use the integration by parts for the integral 
J v^e~^dv = — f ^ dv and obtain the recursive formula 

Iq{a) = a’^~^Ii{a) + {q - 2)/q_2(a) 

that yields the conclusion. □ 

We just present a sketch of the proof of Theorem |2.13[ The detailed computations will be proposed in 
Appendix. The ideas are similar to the proof of Proposition 2.8 and Proposition 2.9 but even more technical 
and laborious. 


Proof of Theorem 2.13 In subsection 2.1 we saw that the transition density of the (/3i,/32)-SBM with drift y 
has an integral representation as in equation [T^ Lemm a |2.14| gives us the expression of the Green functions. 
One can make the change of variable (j){X) =\/2\ + y? proceeding as in Figure [^a. We can show that zero 
is an erasable singularity for the integrand, that is also holomorphic on the entire imaginary line. Since we 
assumed Piy > 0, P 2 y- > 0, the integrand has no poles in (0,^^]. Therefore, being in the case of Figure [^bl 
and since an analogous of Lemma |2.7| holds, one can deform the contour to the imaginary line. One obtains the 
transition density as 


pi^^’^^Ht,x,y) = -e "“s 


— _p-‘^t+v(y-x) 


274 






PiP 2 e-’^™^{w‘^ + y?) + {w - iPiy){w - iP 2 y) 


dw. 
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If 7/; 0 then 




< 1 hence the transition density can be written as 


„(/3i,/32) 

Ffi 




e 2 


27T 


“'R k=o “ ^PlN{w - iP2P)\ 


where we can exchange integral and limit because the series of the absolute values is integrable. We now interpret 
the expression for as a series of Fourier transforms 




Fj^k '■= F 1-^ e ”2 Cj{y, y\/i; iw){w‘^ + yPtY 


The Fourier transform Fj^ki^j,k) can be rewritten as the convolution of Fourier transforms 


0 


Fj,kiuJj,k) = (e ^ 


-1 


{w — iPiyL\/i)^+^{w — iP 2 y'/i)^~^^ 


( 22 ) 


(23) 


-1 


The Fourier transform F , 7 --=--t-t^x 

using jointly Lemma |2.15| and Lemma |2.16[ One concludes the proof using the 


is computed using Lemma 


2.16 


if Pi = P 2 , otherwise 
properties of the i terated 


derivatives of Gaussian densities, and introducing Jq{uj,A) = e 2 +"^‘^J^(—(oj -(- (see Lemma 
more details see the Appendix. 


2.171) • For 

□ 


We assumed both Piy and P 2 fJ. to be positive because, if Piy < 0 or /32/r < 0, the exact computation for 
the density can be even more subtle. This is due to the possible presence of an additional term in the contour 
integral corresponding to zeros of the denominator lying the positive real semi-axis (as in the case of the /3i-SBM 
with drift /r, see Figure [^b2). These cases will be treated in a incoming paper on the exact simulation of a 
Brownian diffusion with drift with several discontinuities [1]. 

Another possible approach in order to solve the § could be to apply the technique used in in case of 
Brownian motion with drift between two barriers, but our approach seems to be more fruitful. 


2.4 Limit cases 


For particular choices of the parameters formulas (20) and (21) reduce to the more simple cases studied before. 

For P 2 — 0, the correspondent barrier Z 2 is completely permeable, so it is like if it disappears, hence one 
would expect to obtain the density of the /3i-SBM with drift. 

Without directly substituting /32 = 0 in the final expression of the transition density, one can notice in equation 


(22) that only Fj^o{ujjfi) for j G {1,2} do not vanish. Moreover equation (22) turns out to be equation (19) with 


P = Pi such that Piy > 0. 

Even for Z 2 —>■ +00 one would expect to obtain the density of the /li-SBM with drift. In fact if the second 
barrier is very far from the starting point of the process, at every finite time the trajectory has no way to see 
the latter barrier and is effected only by the reflection coefficient Pi. 

Less heuristical and more direct would the following approach. First notice that, since Z 2 —>■ -boo, a 3 {x, y), 04(0;, y) 
and z go to -l-oo which implies ujj^k — >■ 00 as soon as A: 7b 0 or j 7b (1^ 2}. Then consider the expression for Fj^k'- 
in equation (22), it is a Fourier transform of a L^-function, hence it is in Lp. It can be shown that it admits 


a limit at infinity, hence this limit has to be zero. Therefore the not vanishing terms in equation are again 
given by j = 1,2 and k = 0. 


3 Exact simulation 

To simulate exactly a process means to simulate it from its law sampling exactly from its finite dimensional 
distributions without approximations (beyond the machine’s). Exact sampling of a random variable can be 
achieved using the rejection sampling method, introduced in |22j . 

The rejection method allows to sample from the density h of a random variable X ^ h{x)dx knowing how 
to sample another one Y ~ g{x)dx if h < Mg ioi M a finite strictly positive constant. The sample y = Y 
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is accepted as a sample of X if and only if m < 


Hv) 

Mg{y) 


where u is the sample of a uniform random variable 

hjy) 


U ~ ^ 0 , 1 ]- Notice that >‘(v) j is a Bernoulli random variable with random parameter Moreover 

the densities h{x) and g{x) do not need to be normalized. 

In our framework, the one-dimensional projection at time < of a (/3i, / 32 )-SBM has a density whose ratio with 
respect to the well known transition probability density of the Brownian motion is a series, as we saw in equation 
(15). What happens if the density cannot be evaluated exactly, since it is an infinite sum? The technique we 
are going to propose allows to evaluate only a finite number of terms of the series, and at the same time, to 
maintain the exactness of the sampling. 


3.1 Generalized rejection sampling method 

Let us introduce our method by explaining a toy example for simulating exactly a Bernoulli random variable 

X Bp with unknown parameter p G [0,1]. If the parameter is known, then clearly X hence an 

exact simulation consists in sampling the uniform random variable U ^ ^[o,i] checking if the sample is 
smaller (or bigger) than p to decide if X = 1 (or X = 0 ). 

Lemma 3.1. Suppose p is an unknown parameter which is approximated by a sequence {pn)n and the rate of 
convergence is at least (Sn)n where (^„)„ is a decreasing vanishing sequence (i.e. \p — Pn\ < Sn)- Then it is 
possible to simulate exactly a Bernoulli of parameter p since X := 1{3„; |f 7 -p„|>i 5 „. c/<p„} ~ I3p . 

Proof. First of all we need to show that, a.s., there exists an n such that \U — Pn\ > dn- Notice that a.s. 
\U — p\ > 0. Since 5n —t 0, a.s. there exist no such that \U — p\ > 25no- Therefore there exist an n (for 
monotonicity it works for n > no) such that a.s. \U — p„| > (5„. 

Now, since 

{U <p} = {3n e N; \U - p„| > 6n, U < p„} = [J {U < pn - Sn} , 

nGN 

then p = P([/ < p) = P(X = 1). 


0 


-V- 

U -|- Sn 

Pn 

£ -u' - (5„ 


1 




1- * - 

Pn ^n 


Pn ^n 



1 U P 

-^-•- 

U + 

- •- 

Pn 

-#- 

1 

0 

Pn 

-^-#- 

u' - Sn 
-•- 

P U 1 ^ 

Pn 



-1 

Pn H" ^n 


1- 

Pn ^n 


— — 1 ^ 1 

Pn H" ^n 


Figure 3: The pictures illustrate the way to sample a Bernoulli random variable X of unknown parameter p\ if 
u < p then X := 1, otherwise X := 0. In the first image u < Pn — Sn hence u < p (resp. u' > Pn + Sn hence 
u' > p). The second images show that, if u < Pn — Sn < u + Sn < Pn (resp. Pn < u' — Sn < Pn + Sn < u) then 
u < p (resp. u' > p) anyway. 


The scheme of the algorithm then will be: 

1 . sample from lA^ we obtain u 

2 . find n such that \u — p„| > Sn, 

3. if M < Pn, then u <p hence X := 1 otherwise X := 0 (see Figure [^. 

□ 

This idea allows us to extend the rejection sampling method for sampling X ^ h{x)dx knowing an approxi¬ 
mation of the density h{x). 
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Theorem 3.2 (Generalized rejection sampling method). Assume one knows how to sample the random variable 
Y with (unnormalized) density g(x). Then one can sample the random variable X with (unnormalized) density 
h{x) under the following assumptions: 

(i) the ratio between the functions g and h is bounded: 


3M > 0 such that 0 < f{y) := 


1 Hy) 

M g{y) 


< 1 for all y gM.] 


(ii) there exists a sequence of explicitly computable functions {fn)n converging to f at a decreasing explicitly 
computable rate (Sn)n- 


Then X ^ (F| 3n; U < fniX) — 6n) i-e. an exact simulation is possible. 


Proof. It is well known from the standard rejection sampling that X ^ {Y\U < f{Y)) (see for example |2D])- 
Lemma 3.1 ensures that we can simulate exactly without knowing f{Y) with complete accuracy. The accept¬ 
ability of the draw y = Y as a sample from X is a Bernoulli with parameter /(y) and we can compute explicitly 
a sequence converging to this quantity fn(y){='- Pn) and its rate of convergence {Sniy))n- Thus the rejection 
sampling scheme based on Lemma |3.1| is the following 


1. sample u from a uniform random variable U ~ ^o,i]j 

2. sample from the density g: we get y = Y, 

3. take y as a sample of X if u < f{y), otherwise reject and start again. More precisely: 

3a. find n such that |w — /n(y)| > 

3b. check whether /„(y) > u, 

3c. if yes accept X = y, if not reject it. 


□ 


3.2 Sampling from the density of the (/3i,/? 2 )-SBM 


We now apply Theorem 3.2 for sampling from the density at time t of the (/3i,/32)-SBM starting at x. We 
already noticed in Proposition |2.8| that its density is absolutely continuous with respect to the one of the 
Brownian motion (t, x, y) with ratio 


v^^^’^^Ht,x,y) = ^(-/3i/32)''^Cj(y,/3i,/32)e’ 

i=i 


{a-j (x,y)+2zk) 




(x,y)+2zk 


(24) 




Since we can only evaluate the sum of the series v^^^’^‘^\t,x,y) with some error, we check if the hypothesis (i) 
and (ii) of Theorem 3.2 are satisfied. 

Lemma 3.3. There exists an upper bound for x,y) uniform in x and y: 


sup 

x,yGM. 


a;, y) 


< V := 


(l + l/?l|)(l + l/?2|) 

1-|/3i/32| 


Proof. 




E'=iIg(2/>^i./32)I (1 + |/3i|)(1 + |/32|) 




□ 
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We denote the truncated series at the first N terms by 

N 4 

:= y](-/3i/32)''y]c4y,/3i,^2)e“ 

fc—0 i—1 


(aj{x,y)-\-2zk)‘^ a j {x ,y)-\-2zk 

— -2t--i-, 


and the rest by {t, x,y) = v^^^’^^\t,x,y) — v]^^’^^\t,x,y). 

Lemma 3.4. The rest of the truncated series is bounded uniformly in x,y: 

Proof. 

< ( y] |/?i/ 32|'') [y]|c4y,/3i,/32)| ) =iJ |^i/32| 


N+l 


\k=N+l 


0 = 1 


□ 


We can apply Theorem |3.2| with 


fn{y) ■= - 


= - 




{t,x,y) and <5„ := |^i^2| 


n+l 


We are then able to sample from the density y i—)■ p^^^'^^\t,x,y) in equation (15) through the generalized 
rejection sampling algorithm and therefore we are able to simulate exactly the Markov process (/3i+2)-SBM 
(for example see Figure]^. 

To increase the efficiency of the rejection algorithm, we apply the following principle: assume we have just 


computed 


0,t,x 


J N 


{y)-' 


and noticed that it is smaller than +, we then take the first index N greater than 


the quantity (log+) ^ log 
algorithm in case it does not find the desired conditions 3a. in Theorem 3.2 


Moreover it is better to fix an integer Nmax in order to stop the 

This index Nmax should be such 


that the r est o f the series is sufficiently small for considering the truncated sum as a good approximation (due 
to Lemma 3.4 an upper bound for the error is v |/3i/32+™“”’). In any case the simulation turns out to be always 
exact (that is the acceptance or rejection is obtained for an index smaller than Nmax) if l/^i+l is not too close 
to 1. In that case we may increase the index Nmax in such a way that is small. 




Figure 4: Exact simulation of a path of the 
(0.7, —0.2)-SBM starting at time 0 in cc = —0.3. 
The barriers are zi = 0 and Z2 = I- 


Figure 5: Comparison between the function y !—>■ 
j)(2--2)(l,0.5,y) obtained from 50000 exact sim¬ 
ulations through generalized rejection sampling 
method and the plot of its truncated version at 
the tenth term {Nmax = 10). The barriers are 
Zi = 0 and Z 2 = 1. 
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Let us compare now the approximation of the density y i— 2 ;, J/) in equation (151 obtained truncating 
the series at the N^ax-^^ term and an histogram of a large number of exact samples from the untruncated density 
computed through the generalized rejection sampling method. For simplicity, we always take time t = 1, starting 
point X = 0.5 and we assume that the barriers are fixed in = 0 and Z 2 = 1- 

We represent in Figureas typical situation, the function y >-)• (1^ 0.5, y). In this case, 100% of the 

50000 simulations are exact. The average number of terms of the series that are necessary in order to decide if 
to accept or reject the simulations is smaller than 2 (1.6). From now on we will denote this number as Nj-ej- 
The transition density in this case is mainly concentrated inside of the interval between the barriers {zi,Z 2 ) 
since /3i > 0 and f32 < 0. Choosing Nmax = 10 the truncated series differs from the untruncated one at most of 
v\j3ij32\^^ ~ 6 • 10 “’’. 




Figure 6 : Comparison between the function y !—>■ 
p(o.3.-o.7)^2,0.5, 2 /) obtained from 50000 exact sim¬ 
ulations through generalized rejection sampling 
method and the plot of its truncated version at 
the tenth term {N^ax = 10). (5io = 3.5 • 10“® and 
Bv = 2.8. The barriers are Zi = 0 and Z 2 = 1. 
The average acceptance number is Nrej = 1.28. 


Figure 7: Comparison between the function y 1 —> 
p(-o. 7 ,o. 3 ) obtained from 50000 exact sim¬ 

ulations through generalized rejection sampling 
method and the plot of its truncated version at 
Nmax — 10. The barriers are zi = 0 and Z 2 = 1. 
The average acceptance number Nrej is 1.27. 


In Figure and we propose skewness parameters with different absolute values and pointing respectively 
inward and outward. All our simulations are exact and Nrej ~ 1-3 is low as expected. In these cases Sn = 0.21""''’ 
and V = 2.8. We can observe in Figure that the process tends to stay between the barriers because when it 
reaches the barrier zi it has probability —^ = 0.65 to be reflected to this region and when it reaches Z 2 the 
probability is 1-^ = 0.85. If the process leaves {zi,Z 2 ), then the probability to be before Zi is larger than to 
be after Z 2 because 1 — /3i > 1 -I- /32. 

In Figure the parameters /3i = —0.7 and /32 = 0.3 induce that the process is more likely to be outside the 
region between the barriers because it is reflected outside this region with probability = 0.85 in zi and 
with probability 0.65 in Z 2 . 
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Density from exact simulation histogram 
D itsinterpolation 

• Approximation truncated at N_{max} 


Density from exact simulation histogram 
D itsinterpolation 

• Approximation truncated at N_{max) 



Figure 8: Comparison between the function y !—>■ 
p(-o.8,-o.6)^2,0.5, 2 /) obtained from 50000 exact 
simulations through generalized rejection sampling 
method and the plot of its truncated version at 
Nmax = 10 with (5io = 3.12 • lO”'^ (and v = 5.54). 
The barriers are zi = 0 and Z 2 = 1. The average 
acceptance number is 3.58. 


0.9 -■ 



Figure 9: Comparison between the function y !—>■ 
0.5,1/) obtained from 50000 exact sim¬ 
ulations through generalized rejection sampling 
method and the plot of its truncated version at 
^max = 20. The barriers are Zi = 0 (completely 
reflecting) and Z2 = 1 (semipermeable). The aver¬ 
age acceptance number is 2.36. 


Figure 1^ represents a case of /3i/32 > 0. From the simulated density function it is confirmed the behaviour 
we expected: the process after a time t will be more likely to stay on the the left (respectively right if the 
parameters are positive) side of the barriers. We chose the parameters Pi < P 2 in such a way that the process 
would more likely stay in (—00, Zi). 

Another interesting example is the case of a completely reflecting barrier and a partially reflecting one: in 
Figure]^ we choose /3i = 1 and P 2 < 0, i.e. zi totally reflecting and Z 2 semipermeable with semipermeabiliy 
coefficient P 2 = —0.4. The process shows the tendency to stay in the between the barriers (zi, Z 2 ), while it will 
have probability zero to be in (—00, Zi). 


Appendix: details in the proof of Theorem 2.13 


We now propose with more details the steps between the convolution of Fourier transforms (221 and the final 
result of Theorem 12.131 

Equation (23) is the convolution of two Fourier transforms, hence one needs to compute first the Fourier 
transforms separately and then the convolution. 

2 

Lemma 3.5. The Fourier transform T of the function w ^ e~^Cj{y, iw){w‘^ + yft)^ is 


h—0 m—0 ^ ^ 


^m+h 

^y2m-\-h 


e " , 


where the functions cj and cj^h for j = 1, 2,3,4, h = 0,1, 2 are given in Lemma 2.14 
Proof. Simply recall that 

2 

h—{) 


Cj{y,yVt;iw) = ^Cj^ 2 -h{y)i^w^ 


and that 


m=0 


(25) 
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Finally one computes the Fourier transforms 


j 2 m-\-h 2 ^ 2 m+h 2 




^y2m-\-h 


e 2 


□ 


(26) 


and concludes. 

If Pi ^ P 2 , as corollary of Lemma 2.15 and Lemma 2.16 one has 
J- (w !—>■ -— I (w) = 

V {w — i/ 3 i/x-\/t)^+i(w — J 

^ ^ - A)"(.v^)"-" [9 (-,amv^) - (-i)-9(-. ft.v^) 

where we defined the function 

g{u;, A) := (21r+(A) — 1) ((21r+(^) — 1) ui) 

Since /3i/r, are both positive then g{u), Pigi/t) = (w), but we give here the proof in the general 

case PigLi/i ^ 0. 


If Pi = P 2 Lemma 
Let us define 


2.16 


gives the formula for T ■ 


/ j 2 m-\-h 2 \ 

■= (-1)'^ * ^^2m+h ^~'^) H 

where A is a fixed real parameter and also 


Fj^k in equation (23) is given by 

p. , _ (-l)’”(2fc-n)! _ y^ 2 r,,,/7i2-h -c o U 

— l^n =0 l^m =0 (k-n)\(k-m)\n\m\k\ C/?, ^vd-/32A«vd)2‘^+i-" ^h=o''j (/j n pi 7= p2, 

if /3l = P2- 


-0 Z.^ITI=0 (fc-n)!(fc—m)!n!7Ti!fc! -, * - ^ 


m,2fc+l ’ 


It remains to compute the function ^{oj, A). One can use that 


fjn 2 2 

where Hn{w) are the Hermite polynomials. 

* ■f^2m+;i(w)e"'^) (w) 

= (21r+(A) - 1 ) / Ir-( v 4 w)?w”e“ ' E +^'"iL2m+?i(w - w)dw 

Jr 

{-w-(uj + A))^ 

lR-(Aw)w”e 2 H2m+h{^ - w)dw 

O M 

(21r+ (yl) — 1) e 2 f (^A(v + A + ^)){v + A + uPf^e 2 H2m+h{~ A — v)dv 

Jr 

= (21r+(A) — 1) e 2 f {v + A + u))'^e ^ H2m+h{~ A — v)dv 

•^(-(21,8+-1)0° -(‘^+^1) 

We can then use the binomial formula and then the explicit expression for the Hermite polynomials: 

LtJ 


= (2lR+(H)-l)eE+^- 
{y—w — A—(j 




w 


. 71-21 


21 





















Therefore 


— (21 r+ (^) — 1 ) e 2 
2 +LIJ 


-^+Au; 


(v + A + uj)^e 2 H 2 m+h(~(A -h v))di 


= (2m + h)! g £!(2(m - £) + h)! 

where, if /(w, A) denotes the interval (— (21]i{+(yl) — 1) oo, — (w + A)), 

= (-l)'^(2lK+(A)-l)e"^+^“ [ {v + A + uj)^e-'^{A + vf^"^-^'>+'^dv 

J lioj^A) 


n 2{m—i)+h 


r=0 s=0 \ ® / 

where = (21 r+(A) - 1) = 


^+^‘"/,+,(-(a; + A)) A >0 . 


-e^+^-4+«(-(w + A)) A<0 


with /(j 


and Iq defined in Lemma 
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